"""
问题一简化分析脚本
只保存图片，不显示，避免交互问题
"""

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import os

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

def load_and_analyze_spd():
    """加载并分析SPD数据"""
    try:
        # 从CSV文件读取
        df = pd.read_csv('Problem 1.csv', encoding='gbk')
        wavelengths = []
        spd = []
        import re
        for i, row in df.iterrows():
            wl_match = re.search(r'(\d+)', str(row.iloc[0]))
            if wl_match:
                wl_value = float(wl_match.group(1))
                intensity_value = float(row.iloc[1])
                if 380 <= wl_value <= 780:
                    wavelengths.append(wl_value)
                    spd.append(intensity_value)
        wavelengths = np.array(wavelengths)
        spd = np.array(spd)
        print("✓ 成功从CSV文件读取数据")

        return wavelengths, spd

    except Exception as e:
        print(f"✗ 数据读取失败: {e}")
        return None, None

def calculate_spectral_parameters(wavelengths, spd):
    """计算光谱参数"""
    import colour

    # 构造光谱对象
    spd_colour = colour.SpectralDistribution(dict(zip(wavelengths, spd)), name='LED SPD')

    # 计算XYZ三刺激值
    XYZ = colour.sd_to_XYZ(spd_colour)

    # 计算色品坐标
    xy = colour.XYZ_to_xy(XYZ)

    # 计算相关色温
    CCT = colour.xy_to_CCT(xy, method='McCamy 1992')

    # 计算Duv
    uv = colour.XYZ_to_UCS(XYZ)
    u = (4 * uv[0]) / (uv[0] + 15 * uv[1] + 3 * uv[2])
    v = (6 * uv[1]) / (uv[0] + 15 * uv[1] + 3 * uv[2])

    from colour.temperature import CCT_to_uv
    uv_bb = CCT_to_uv((CCT, 0))
    Duv = np.sqrt((u - uv_bb[0])**2 + (v - uv_bb[1])**2)

    # 计算CIE Ra
    from colour.quality import colour_rendering_index
    Ra = colour_rendering_index(spd_colour)

    # 计算mel-DER
    from colour import MSDS_CMFS
    cmfs = MSDS_CMFS['CIE 1931 2 Degree Standard Observer']
    photopic_wavelengths = cmfs.wavelengths
    photopic_sd = cmfs.values[:, 1]
    photopic_interp = np.interp(wavelengths, photopic_wavelengths, photopic_sd)

    # 读取melanopic数据
    melanopic_wavelengths = []
    melanopic_values = []
    try:
        with open('CIE S 026E.csv', encoding='gbk') as f:
            import csv
            reader = csv.reader(f)
            next(reader)
            next(reader)
            for row in reader:
                if len(row) < 2:
                    continue
                wl = row[0].strip().replace('"','').replace(' ','')
                val = row[1].strip().replace('"','').replace(' ','').replace('E','e').replace(',','.')
                try:
                    wl = float(wl)
                    val = float(val)
                    melanopic_wavelengths.append(wl)
                    melanopic_values.append(val)
                except:
                    continue
    except:
        print("⚠ 无法读取melanopic数据")

    if len(melanopic_wavelengths) > 0:
        melanopic_wavelengths = np.array(melanopic_wavelengths)
        melanopic_values = np.array(melanopic_values)
        melanopic_interp = np.interp(wavelengths, melanopic_wavelengths, melanopic_values)
        mel_DER = np.sum(spd * melanopic_interp) / np.sum(spd * photopic_interp)
    else:
        mel_DER = 1.0

    return {
        'XYZ': XYZ,
        'xy': xy,
        'CCT': CCT,
        'Duv': Duv,
        'Ra': Ra,
        'mel_DER': mel_DER,
        'photopic_interp': photopic_interp,
        'melanopic_interp': melanopic_interp if len(melanopic_wavelengths) > 0 else None
    }

def create_separate_visualizations(wavelengths, spd, params):
    """创建分别的可视化图表"""

    print("正在生成可视化图表...")

    # 1. 光谱分布与响应函数对比
    plt.figure(figsize=(12, 8))
    ax1 = plt.gca()
    ax1_twin = ax1.twinx()

    ax1.plot(wavelengths, spd, 'b-', linewidth=2, label='LED光谱', alpha=0.8)
    ax1.fill_between(wavelengths, spd, alpha=0.3, color='blue')
    ax1_twin.plot(wavelengths, params['photopic_interp'], 'g-', linewidth=2,
                  label='V(λ)光视效率', alpha=0.8)
    if params['melanopic_interp'] is not None:
        ax1_twin.plot(wavelengths, params['melanopic_interp'], 'r-', linewidth=2,
                      label='黑视响应', alpha=0.8)

    ax1.set_xlabel('波长 (nm)')
    ax1.set_ylabel('光谱功率密度 (W/nm)', color='b')
    ax1_twin.set_ylabel('响应函数', color='g')
    ax1.set_title('光谱分布与响应函数对比')

    lines1, labels1 = ax1.get_legend_handles_labels()
    lines2, labels2 = ax1_twin.get_legend_handles_labels()
    ax1.legend(lines1 + lines2, labels1 + labels2, loc='upper right')

    plt.tight_layout()
    plt.savefig('1_光谱分布与响应函数.png', dpi=300, bbox_inches='tight')
    plt.close()
    print("✓ 1_光谱分布与响应函数.png")

    # 2. 色度图
    plt.figure(figsize=(10, 8))
    ax2 = plt.gca()

    x_boundary = [0.175, 0.175, 0.734, 0.734, 0.175]
    y_boundary = [0.005, 0.816, 0.816, 0.005, 0.005]
    ax2.plot(x_boundary, y_boundary, 'k-', linewidth=1)

    ax2.plot(params['xy'][0], params['xy'][1], 'ro', markersize=12,
             label=f'LED光源\n({params["xy"][0]:.3f}, {params["xy"][1]:.3f})')

    try:
        from colour.temperature import CCT_to_xy
        cct_range = np.linspace(4000, 25000, 50)
        bb_x, bb_y = [], []
        for cct in cct_range:
            try:
                xy_bb = CCT_to_xy(cct)
                bb_x.append(xy_bb[0])
                bb_y.append(xy_bb[1])
            except:
                continue
        if len(bb_x) > 0:
            ax2.plot(bb_x, bb_y, 'k--', linewidth=1, alpha=0.7, label='黑体轨迹')
    except:
        pass

    ax2.set_xlabel('x')
    ax2.set_ylabel('y')
    ax2.set_title('CIE 1931色度图')
    ax2.legend()
    ax2.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.savefig('2_CIE1931色度图.png', dpi=300, bbox_inches='tight')
    plt.close()
    print("2_CIE1931色度图.png")

    # 3. 参数雷达图
    plt.figure(figsize=(10, 8))
    ax3 = plt.subplot(111, projection='polar')

    categories = ['CCT\n(归一化)', 'CIE Ra', 'mel-DER', 'Duv\n(归一化)']

    cct_norm = min(100, max(0, (params['CCT'] - 2000) / (8000 - 2000) * 100))
    duv_norm = max(0, 100 - params['Duv'] * 1000)
    ra_norm = params['Ra']
    mel_norm = params['mel_DER'] * 100

    values = [cct_norm, ra_norm, mel_norm, duv_norm]

    angles = np.linspace(0, 2 * np.pi, len(categories), endpoint=False).tolist()
    values += values[:1]
    angles += angles[:1]

    ax3.plot(angles, values, 'o-', linewidth=2, color='red')
    ax3.fill(angles, values, alpha=0.25, color='red')
    ax3.set_xticks(angles[:-1])
    ax3.set_xticklabels(categories)
    ax3.set_ylim(0, 100)
    ax3.set_title('光源参数雷达图', pad=20)

    plt.tight_layout()
    plt.savefig('3_参数雷达图.png', dpi=300, bbox_inches='tight')
    plt.close()
    print("3_参数雷达图.png")

    # 4. 参数对比柱状图
    plt.figure(figsize=(12, 8))
    ax4 = plt.gca()

    param_names = ['CCT (K)', 'CIE Ra', 'mel-DER', 'Duv']
    param_values = [params['CCT'], params['Ra'], params['mel_DER'], params['Duv']]
    colors = ['#FF6B6B', '#4ECDC4', '#45B7D1', '#96CEB4']

    bars = ax4.bar(param_names, param_values, color=colors, alpha=0.7)
    ax4.set_ylabel('数值')
    ax4.set_title('光源参数对比')
    ax4.tick_params(axis='x', rotation=45)

    for bar, value in zip(bars, param_values):
        height = bar.get_height()
        ax4.text(bar.get_x() + bar.get_width()/2., height + max(param_values)*0.01,
                f'{value:.2f}', ha='center', va='bottom')

    plt.tight_layout()
    plt.savefig('4_参数对比柱状图.png', dpi=300, bbox_inches='tight')
    plt.close()
    print("4_参数对比柱状图.png")

    # 5. 光谱功率分布
    plt.figure(figsize=(12, 8))
    ax5 = plt.gca()
    ax5.plot(wavelengths, spd, 'b-', linewidth=2)
    ax5.fill_between(wavelengths, spd, alpha=0.3, color='blue')
    ax5.set_xlabel('波长 (nm)')
    ax5.set_ylabel('光谱功率密度 (W/nm)')
    ax5.set_title('LED光源光谱功率分布')
    ax5.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.savefig('5_光谱功率分布.png', dpi=300, bbox_inches='tight')
    plt.close()
    print("5_光谱功率分布.png")

    # 6. 色温评价
    plt.figure(figsize=(10, 8))
    ax6 = plt.gca()

    cct_categories = ['暖白光\n(<3000K)', '中性白\n(3000-5000K)', '冷白光\n(>5000K)']
    cct_values = [0, 0, 0]

    if params['CCT'] < 3000:
        cct_values[0] = 1
    elif params['CCT'] < 5000:
        cct_values[1] = 1
    else:
        cct_values[2] = 1

    colors_cct = ['#FFB347', '#87CEEB', '#E6E6FA']
    bars_cct = ax6.bar(cct_categories, cct_values, color=colors_cct, alpha=0.7)
    ax6.set_ylabel('分类')
    ax6.set_title(f'色温分类 (CCT={params["CCT"]:.0f}K)')
    ax6.set_ylim(0, 1.2)

    plt.tight_layout()
    plt.savefig('6_色温评价.png', dpi=300, bbox_inches='tight')
    plt.close()
    print("6_色温评价.png")

    # 7. 显色性评价
    plt.figure(figsize=(10, 8))
    ax7 = plt.gca()

    ra_categories = ['较差\n(<60)', '良好\n(60-80)', '优秀\n(≥80)']
    ra_values = [0, 0, 0]

    if params['Ra'] < 60:
        ra_values[0] = 1
    elif params['Ra'] < 80:
        ra_values[1] = 1
    else:
        ra_values[2] = 1

    colors_ra = ['#FF6B6B', '#FFD93D', '#6BCF7F']
    bars_ra = ax7.bar(ra_categories, ra_values, color=colors_ra, alpha=0.7)
    ax7.set_ylabel('分类')
    ax7.set_title(f'显色性评价 (Ra={params["Ra"]:.1f})')
    ax7.set_ylim(0, 1.2)

    plt.tight_layout()
    plt.savefig('7_显色性评价.png', dpi=300, bbox_inches='tight')
    plt.close()
    print("7_显色性评价.png")

    # 8. 生理节律影响
    plt.figure(figsize=(10, 8))
    ax8 = plt.gca()

    mel_categories = ['较弱\n(<0.5)', '较强\n(≥0.5)']
    mel_values = [0, 0]

    if params['mel_DER'] < 0.5:
        mel_values[0] = 1
    else:
        mel_values[1] = 1

    colors_mel = ['#98D8C8', '#F7DC6F']
    bars_mel = ax8.bar(mel_categories, mel_values, color=colors_mel, alpha=0.7)
    ax8.set_ylabel('分类')
    ax8.set_title(f'生理节律影响 (mel-DER={params["mel_DER"]:.3f})')
    ax8.set_ylim(0, 1.2)

    plt.tight_layout()
    plt.savefig('8_生理节律影响.png', dpi=300, bbox_inches='tight')
    plt.close()
    print("8_生理节律影响.png")

    # 9. 各波段功率分布
    plt.figure(figsize=(12, 8))
    ax9 = plt.gca()

    blue_power = np.sum(spd[(wavelengths >= 400) & (wavelengths < 500)])
    green_power = np.sum(spd[(wavelengths >= 500) & (wavelengths < 600)])
    red_power = np.sum(spd[(wavelengths >= 600) & (wavelengths < 700)])

    bands = ['蓝光\n(400-500nm)', '绿光\n(500-600nm)', '红光\n(600-700nm)']
    band_powers = [blue_power, green_power, red_power]
    colors_bands = ['#4169E1', '#32CD32', '#DC143C']

    bars_bands = ax9.bar(bands, band_powers, color=colors_bands, alpha=0.7)
    ax9.set_ylabel('功率 (W)')
    ax9.set_title('各波段功率分布')
    ax9.tick_params(axis='x', rotation=45)

    for bar, value in zip(bars_bands, band_powers):
        height = bar.get_height()
        ax9.text(bar.get_x() + bar.get_width()/2., height + max(band_powers)*0.01,
                f'{value:.2f}', ha='center', va='bottom')

    plt.tight_layout()
    plt.savefig('9_各波段功率分布.png', dpi=300, bbox_inches='tight')
    plt.close()
    print("9_各波段功率分布.png")

    # 10. 色品坐标分析
    plt.figure(figsize=(10, 8))
    ax10 = plt.gca()

    ax10.scatter(params['xy'][0], params['xy'][1], s=200, c='red', alpha=0.7)
    ax10.set_xlabel('x')
    ax10.set_ylabel('y')
    ax10.set_title('色品坐标分析')
    ax10.grid(True, alpha=0.3)
    ax10.text(0.05, 0.95, f'x={params["xy"][0]:.4f}\ny={params["xy"][1]:.4f}',
              transform=ax10.transAxes, verticalalignment='top',
              bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))

    plt.tight_layout()
    plt.savefig('10_色品坐标分析.png', dpi=300, bbox_inches='tight')
    plt.close()
    print("10_色品坐标分析.png")

    # 11. XYZ三刺激值
    plt.figure(figsize=(10, 8))
    ax11 = plt.gca()

    xyz_names = ['X', 'Y', 'Z']
    xyz_values = params['XYZ']
    colors_xyz = ['#FF6B6B', '#4ECDC4', '#45B7D1']

    bars_xyz = ax11.bar(xyz_names, xyz_values, color=colors_xyz, alpha=0.7)
    ax11.set_ylabel('三刺激值')
    ax11.set_title('CIE 1931 XYZ三刺激值')

    for bar, value in zip(bars_xyz, xyz_values):
        height = bar.get_height()
        ax11.text(bar.get_x() + bar.get_width()/2., height + max(xyz_values)*0.01,
                f'{value:.1f}', ha='center', va='bottom')

    plt.tight_layout()
    plt.savefig('11_XYZ三刺激值.png', dpi=300, bbox_inches='tight')
    plt.close()
    print("11_XYZ三刺激值.png")

    # 12. 综合评价
    plt.figure(figsize=(10, 8))
    ax12 = plt.subplot(111, projection='polar')

    evaluation_categories = ['色温\n稳定性', '显色性', '生理\n节律', '光谱\n完整性']

    temp_score = max(0, 100 - abs(params['Duv']) * 10000)
    ra_score = params['Ra']
    mel_score = min(100, params['mel_DER'] * 100)
    spectrum_score = min(100, len(wavelengths) / 4)

    eval_values = [temp_score, ra_score, mel_score, spectrum_score]

    angles_eval = np.linspace(0, 2 * np.pi, len(evaluation_categories), endpoint=False).tolist()
    eval_values += eval_values[:1]
    angles_eval += angles_eval[:1]

    ax12.plot(angles_eval, eval_values, 'o-', linewidth=2, color='purple')
    ax12.fill(angles_eval, eval_values, alpha=0.25, color='purple')
    ax12.set_xticks(angles_eval[:-1])
    ax12.set_xticklabels(evaluation_categories)
    ax12.set_ylim(0, 100)
    ax12.set_title('综合评价', pad=20)

    plt.tight_layout()
    plt.savefig('12_综合评价.png', dpi=300, bbox_inches='tight')
    plt.close()
    print("12_综合评价.png")

    print("\n所有图表已分别保存完成！")

def main():
    """主函数"""
    print("开始LED光源问题一分析...")

    # 1. 加载数据
    wavelengths, spd = load_and_analyze_spd()
    if wavelengths is None:
        return

    # 2. 计算参数
    print("正在计算色度学参数...")
    params = calculate_spectral_parameters(wavelengths, spd)

    # 3. 生成可视化
    create_separate_visualizations(wavelengths, spd, params)

if __name__ == "__main__":
    main()
